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Two of the most popular quantum mechanical models of interacting fermions are compared to 
each other and to potentially exact solutions for a pair of contact-interacting fermions trapped in a 
' ID double-well potential, a model of atoms in a quasi-lD optical lattice or electrons of a Hydrogen 

, molecule in a strong magnetic field. An exact few-body Hamiltonian is solved numerically in mo- 

mentum space yielding a highly-correlated eigenspectrum. Additionally, approximate ground-state 
energies are obtained using both density functional theory (DFT) functional and 2-site Hubbard 
models. A ID adiabatic LDA kernel is constructed for use in time-dependent density functional 
. theory (TDDFT), and the resulting excited-state spectrum is compared to the exact and Hubbard 

\Q ' results. DFT is shown to give accurate results for wells with small separations but fails to describe 

localization of opposite spin fermions to different sites. A locally cognizant (LC) density functional 
based on an effective local fermion number would provide a solution to this problem, and an ap- 
proximate treatment presented here compares favorably with the exact and Hubbard results. The 
TDDFT excited-state spectrum is accurate in the small parameter regime with non-adiabatic ef- 
fects accounting for any deviations. As expected, the ground-state Hubbard model outperforms 
DFT at large separations but breaks down at intermediate separations due to improper scaling to 
the united-atom limit. At strong coupling, both Hubbard and TDDFT methods fail to capture the 
appropriate energetics. 
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PACS numbers: 31.15.Ew, 71.15.Mb, 71.10.- 



I. INTRODUCTION 



One of the underlying challenges in computational physics, notably in density functional theory, is the accurate 
and reliable treatment of many-body interactions in matter. The typical interaction-type that most electrons in 
' solids experience is the Coulomb interaction, and this is a pernicious one owing to its long-ranged nature. Many 
\Q , of the difficulties in providing accurate numeric results for Coulomb-interacting systems are the results of the long- 
range. It would be insightful and beneficial to the development of techniques that explicitly describe many-body 
systems if the Coulomb interaction were instead a more local version. One such replacement, the contact interaction 
offers a powerful test case to probe specifically the many-body problem without long-range complications. This is 
particularly important in the analysis and formal improvement of practical methods of electronic structure theory such 
as Hubbard theories and density functional ones. In this paper, we explore a simple model system of fermions that 
interact via a contact interaction. Remarkably, the interacting quantum problem is shown here to reduce to a solvable 
system of integral equations. Additionally, we extend developments in the density functional modeling of contact- 
interacting systems to the time-dependent domain. This work has implications both for the formal development of 
density-functional theory and for the practical numeric analysis of quasi ID systems. 

The contact interaction can arise as a simplified form of 3D interactions in highly confined quasi-lD systems, exper- 
imentally realizable in ID optical lattices [1-3]. In these situations, the contact term is a reasonable approximation to 
the interaction when the ratio of the interaction strength to the transverse well width is large. An additional approx- 
imation made is that the wave- function, although still 3D, factors into transverse and longitudinal components with 
the interaction only affecting the longitudinal part. This decomposition is valid when the transverse and longitudinal 
length scales are drastically different. We consider only the nontrivial longitudinal part. Thus, the 3D interaction 
is replaced by a ID one. Admittedly, a more realistic reduction of the 3D Coulomb problem to an effectively ID 
one is widely debated, and, here, we have chosen a particular form partially motivated by the advantages it offers 
computationally. 

One of the computational advantages of the <5-function interaction within the DFT formalism is that the exchange 
and Hartree density functionals are explicit functionals of the fermion density, a result implying that the local-density 
approximation (LDA) to exchange is self-interaction free. Thus, the functionals used for Hartree and exchange 
represent the exact- exchange formalism (EXX). This locality of EXX is not true for the long-ranged 3D Coulomb 
interaction, and performing EXX in 3D is a significantly more complicated endeavor. For the full LDA in the ID 
case, any inaccuracies are caused by the mismodeling of correlation. In 3D, the inaccuracies are mixed between 
exchange and correlation due to the long-ranged nature of the Coulomb potential. A notable example of a long- 
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ranged correlation problem is the inability of the LDA to localize single electrons on distant sites. Hubbard theory, 
on the other hand, preserves this limit but sacrifices accuracy. 

For this reason and the clarifying simplicity of ID models, the <5-function model for ID fcrmions has been consid- 
ered in several recent studies. In particular, the introduction of a local density correlation functional has provided 
exceptionally accurate results for the ground-state of finite systems such as the Diracium [4] and the ID analog of 
Hooke's atom [4, 5]. A modified parameterization of the local density approximation (LDA) has been used to describe 
interacting fermions in harmonic confining potentials [5]. A time-dependent Thomas- Fermi theory has been used to 
predict the excited-state properties of a many-fermion system but did not amount to the Kohn-Sham (KS) formulation 
of DFT [6] . More general problems involving this contact interaction also have been studied recently [7-9] . 

In this paper, we investigate a two-site <5-function problem, the ID quantum analog of the H2 molecule [10] or a 2-site 
optical lattice, hitherto referred to as ID H 2 . The author introduces a ID contact-interacting time-dependent density 
functional theory (TDDFT) and compares the results to the exact solutions and to the 2-site Hubbard results. In order 
to accurately express the ID adiabatic LDA kernel, it is found necessary to improve the existing parameterization of 
the result for the reference system to include higher-order terms at the low- and high-density limits. The resulting 
functional reproduces the previously known expansion terms reported in Refs. [4, 5] but additionally includes a 
critically important high density term needed to describe the adiabatic LDA kernel. The author develops a numeric 
scheme to obtain the exact bound-state spectrum for ID H2 from a set of ID integral equations. The results are tested 
against a general form of the virial theorem. Numeric integration techniques are developed to handle the oscillatory 
nature of the integrand. In the stretched system, the author contrives a scheme to correct the long-range correlation 
self-interaction error. The results are compared to the well-known two-site Hubbard model that is expected to be 
accurate in the long-separation limit. 

Throughout, we assume that our ID fermions have the same mass as electrons, and we use atomic units (e 2 = h = 
m e = 1) so that all energies are in Hartrees and all lengths in Bohr radii. A realistic correspondence between the 
interaction strength A and the details of a trapping potential is given by Ref. [11]. Here, the external potential and 
interaction strengths arc arbitrarily chosen so that the total energies come out in the typical chemical range for the 
sake of physical intuition. 



II. ONE-DIMENSION CONTACT-INTERACTING FERMIONS IN A DOUBLE- WELL 



In this section, we state the quantum mechanical model of trapped ID fcrmions. From the non-interacting case, we 
anticipate how the many-body wave-functions are structured. 
The ID Hamiltonian for the particle-pair in a ID well is 




(1) 

The terms from left to right are the kinetic term, the contact interaction strength, and the trapping potential with 
a being half the inter- well spacing. Z and A are the relative strength of the local and interaction potentials, i labels 
the fermion number, and Xi is the position of the z-th fermion. The lowest energy eigenfunction or ground-state 
wave-function is assumed to vanish infinitely far from the trapping potential. 

The non-interacting problem has two single-particle eigenvalues: e± = —k\j2 with k± = Z + 
l/(2a)LambertW (±2aZexp(— 2aZ))[8, 12]. The '+' corresponds to the lower energy spatially symmetric (gerade) 
single-particle bonding state and the '— ' corresponds to the anti-symmetric (ungerade) single-particle anti-bonding 
state. The Lambert W(y) function is the principle solution for x of y = xexp(x). The normalized single particle wave- 
functions are <j)±(x) — N± (exp(— k±\x — a\) ± exp(— k±\x + a\ j) where a; is a non-interacting particle?s position, and 

N± equals 1/ ± ^4a + -^j exp(— 2fc±a). While a gerade ground-state exists for all values of aZ, the ungerade 

state only exists for aZ > 1/2. We refer to the large separation case as stretched H 2 and the small separation as 
crushed H 2 tending to the united- atom limit. Table I characterizes the various non-interacting two-particle states and 
defines what is meant by the term designations, So, Si, S2, and T\. 

The virial theorem is important in practical applications where it is often used to verify numeric results. In its 
fundamental form, it states that 2(T) = (xj^V(x)}. The expectation value denotes an average over spatial variables 
of the respective eigen- function solutions of Eq. 1. The modified version, generated through integration by parts and 
typically used for the Coulomb interaction, is 2(T) = — (V(x)) but docs not hold true here due to the fixed location 
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of the external potential. A careful integration by parts shows that a generalized version, 

/oo 
dy r(a,y)d/da^(a,y) 
-OO 

/OO 
dyip*(-a,y)d/daip(-a,y), (2) 
-oo 

is valid for the ground-state of the system. Note that we use y her for the spatial coordinate to avoid confusion with X\ 
and X2- d/daip(a, y) is the average of the right and left derivatives with respect to the first argument. The interaction 
potential part needs no modification. A many-site version of Eq. 2 holds for multiple <5-wells but has more terms on 
the right of the equal sign. The right hand side of the general version vanishes for an isolated 5-well and a periodic 
lattice of wells because in these cases the average derivative is null by symmetry and the form of the virial theorem, 
2(f) = -(V(x)), is justified. 

In the interacting problem, there are only two relevant parameters, A and a. The external potential strength can 
be scaled to unity leaving the energy in units of Z 2 times the atomic unit and length in units of 1/Z times the atomic 
unit. This leaves four regimes to consider. A small and a large is the weakly-interacting two well solution. A large 
and a small is the highly correlated double well solution. For sufficiently small a, the wells merge to the united-atom 
limit. When A and a are large, we have a system where long-range correlations can be important. Finally, A small 
and a small is a regime where the interaction can be treated perturbatively, and the well is almost a single well. Large 
a is the domain of validity of the Hubbard model. 

This paper examines 2 of these regimes in detail: the small width double- well, a = 1, and the large range hopping 
scale, a = 2, both with Z = 1. We note that in the former range, the separation between the wells is large enough to 
exceed the united atom limit yet still preserves the existence of at least two bound states. The latter case (aZ = 2) 
is the regime where Hubbard theory is designed to be maximally valid. For the remainder of the paper, we assume 
that Z = 1. 

The two-site Hubbard will be compared to the more general density functional and exact solutions. The Hubbard 
model itself is important for the conceptual picture it provides and its role in the popular LDA+U method [24, 25]. 
Furthermore, the Hubbard picture is designed to reproduce localization in the stretched limit. A property that is not 
matched in the LDA version of DFT. In the Hubbard approach, the Hamiltonian Eq. 1 is simplified to a Hubbard 
Hamiltonian [26]. The simplification is valid when the wave- function overlap between sites is small as is the case for 
a = 2 but not a = 1. The approximate Hamiltonian can be diagonalized exactly and relies on two parameters, t (the 
hopping term) and U (the on-site repulsion). 

To find t, the hopping term, we consider linear combinations of the non-interacting single particle orbitals that give 
left and right localized fermions: <j)R(x) = {4>+{x) + </)_(ir))/v / 2 and 4>l(x) — (<j)+(x) — 0_ (x))/\/2. The hopping term 
is the projection, on one localized wave- function, of the kinetic energy operator acting on the other localized fermion 
wave-function: t = 1/2 J dx (f> R (x)\7 2 cj) L (x) = l/4(fc"j_ — fc 2 .). Note that we defined t to be positive. This definition 
is chosen to reproduce the proper long range limit. Traditionally, the t term is fixed using non-orthogonal localized 
solutions. The difference is negligible in the large separation limit. For the small separation limit, our definition of t 
gives a different result than the traditional hopping term. This is because the 4>- state is no longer bound for small a, 
and consequently, the maximally localized solutions are not really localized to any one site. Curiously, the convention 
used here more accurately describes the non-interacting case at small separations than the traditional definition. 

The Hubbard Hamiltonian can be exactly diagonalized and has eigenvalues Es = 2e + ^ — j^flQt 2 + U 2 , Es 1 = 
2e + U, and E$ 2 = 2e + ^ + \\f\M 2 ~~\-U 2 ■ Physically the first excited singlet state at large separation, a, corresponds 
to 2 fermions localized on one site. This is the Diracium system and a distant empty site. If we make a correspondence 
between the Hubbard energy of S\ and the total energy of Diracium, U can be expressed as a function of A. The 
total energy for Si is written E = 2e + U with e < 0. The 2-particle allowable range for U from Diracium is to — e. 
Above this limit, at X cr it., U^ it = — e, and the single well no longer binds two fermions. We match the Hubbard U 
value to give exact results given in Refs. [13] and [4]. For example, a U of 0.354 gives the correct energy for Diracium 
with A = 1 and Z = 1. For numerical convenience, we parameterize the interaction energy versus A when Z = 1: 
U(X) w 0.500A - 0.163A 2 + 0.017A 3 + C(A 4 ), valid for A < \ cnt with a max error of about 0.5 %. 

Despite scaling problems, the Hubbard model is well trusted at large a because its computational convenience and 
its facility with handling long-ranged ground-state correlations. LDA does not handle these correlations, and even 
the exact solution method must be carefully formulated in this limit. 
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III. TIME-DEPENDENT DENSITY FUNCTIONAL THEORY 

In this section we generalize work, done by us and others, in ground-state density-functional theory to the time- 
dependent case. This extension requires several improvements upon the ground-state theory. There are four main 
results of this section: 1. the discovery of the importance in TDDFT of the third term in the high density correlation 
energy, 2. the creation of an improved correlation functional that properly describes this term, 3. the proposal of a 
local effective fermion measure based on KS orbitals, and 4. the introduction of a new correlation functional that can 
properly describe spatially separated systems. The results can be found in equations 3, 5, 14, and 15 respectively. 

A general approach to find the ground-state energy and fermion density of ID H2 is to use Density Functional 
Theory (DFT). In ground-state DFT, the details of the external potential are kept, but the many-body interaction 
is transformed to an effective local potential derived from the exchange-correlation density functional [15, 16]. Given 
the exact exchange-correlation functional, DFT would return the exact results for the total energy and density. In 
practice, this exchange-correlation contribution is approximated. An active area of research in the 3D case is to 
improve the accuracy and reliability of approximations to the exchange-correlation functional. 

According to spin-density functional theory [15], the ground-state total energy is a functional of the particle density 
and the local magnetization. In this work, we make the assumption that an axis of magnetization is chosen, and 
the local magnetization is given by ((x) = (rif(x) — ni(x))/(n^(x) + ni(x)) where n-\ and rij arc up- and down-spin 
densities projected on the magnetization axis. The total ground state energy can then be decomposed as follows: 
E[n, (} — T s [n, (] + E H [n] + E xc [n, (] + J dx v cxt (x)n(x) in ID where £ H [n] is the exactly known Hartree or classical 
density-density interaction contribution, v ex t(x) is the given inhomogeneous potential, T s [n,Q is the exactly known 
kinetic energy of non-interacting fermions at a given density, E xc [n, £] is the exchange-correlation energy. The solution 
for the problem of interest is found by studying the Kohn-Sham (KS) system, the non-interacting counterpart to the 
physical system [16]. The spin-densities are obtained from the occupied KS orbitals, n cr (x) = J2i occ I^.^C^)! 2 - 

In the contact ID case, the Hartree and exchange terms, E HX [n^,ni] = A J dx n^(x)ni(x), are known exactly, and 
only the correlation energy functional must be approximated in practice. It is important to note that the Hartree and 
exchange functional for this contact interaction is self-interaction free, so in essence, the pure density functional already 
includes explicitly exact-exchange (EXX). Many of the known practical limitations of 3D DFT can be addressed by 
applying EXX but in 3D, this is computationally demanding. The LDA functionals in this paper include the EXX 
formalism, so the analysis and results here will be useful in the next development stage of DFT where EXX-compatiblc 
correlation is considered in detail. In particular, since the interaction here is local, certain difficult-to-modcl long- 
range aspects of correlation will be isolated. For example, contact-interacting fermions will still exhibit spin-density 
waves and long-ranged entanglement, two problems extremely difficult to model using traditional density functional 
methods. 

The correlation energy functional is often modeled using the correlation energy of a solvable reference system of 
fermions. For the chosen ID interaction, this is the Gaudin-Yang model solved exactly via the Bethe-Ansatz technique 
[17-19]. The correlation energy per particle of the uniform system, e" m ^', can then be parameterized as in Refs. [4] 
and [5] to reproduce the exact curve. The high-density expansion is 

£c(n) = 

-A 2 /24 +A 3 C(3)/(27r 4 n) - A 4 0.00094/n 2 + 0(\ 5 /n 3 ). 

(3) 

C(3) is the Ricmann zeta function evaluated at 3. See appendix A for the calculation of the second term. 
The low-density correlation energy is 

e c (n) = 

-An/4 +n 2 7r 2 /8 - nV log(2)/(3A) + 0(n 4 /\ 2 ). 

(4) 

The low density limit can be understood by noting that the interaction is so strong that it mimics Fermi repulsion, 
so it must cancel the Hartree and exchange terms and add a kinetic-like contribution to the energy. The third term 
is referred to in Ref. [7, 20]. 

A modified parameterization of the correlation energy per particle as a (4,4) Pade is 

umf An 3 + Bn 2 + Cn 

c { ' Dn 3 + En 2 + Fn + 1 

(5) 
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with A = -7.031 951, B = -2.169 922, C = -0.25, D = 168.766 814, E = 77.069 721, and F = 13.614 491. This 
parameterization gives 3 terms in both the high- and low-density expansions of the correlation energy. The error is less 
than 0.5% error in the correlation energy per particle for all densities. We note the remarkable fact that parameters, 
A-F, are determined exactly with no approximate numerical fit. Details are given in the appendix B. 

The LDA correlation energy functional is an integral over local contributions of the reference system's correlation 
energy per particle times the fermions per unit volume, 

E% DA [n,C] = jdx n(x) e« n * (n(x))f(((x)). (6) 

The exact high-density limit for /(£) can be obtained via diagrammatic perturbation theory (Appendix C) and is 
approximately /(C) (1 — C 2 ) used here for all densities. The improved (4,4) Pade parameterization of the LDA 
correlation energy functional is used because the ones given in Ref. [4] and [5] do not accurately reproduce the high 
density correlation kernel needed for TDDFT as we will explain later. 

The DFT solution is obtained through a self-consistent solution of the Kohn-Sham equations using a modified 
version of the DFT code in Ref. [4]. In this code, a Numcrov integration scheme is combined with the shooting 
method to obtain solutions of the Kohn-Sham equations. The number of grid points is chosen to converge energies to 
within mHartree accuracy, and the output is checked against an analytic EXX solution given in appendix D. 

Time-dependent DFT allows the determination of the excited states. The exact excitations occur at the poles 
of the density response- function [21, 22]. Finding the excited-state transition energies amounts to the solution of a 
generalized eigenvalue problem: 

ftij<T,klT(V(I))F(I),klT = ^\l) F {I),ija ( 7 ) 

with 

QijaMr(u) = S aT S lk S : ji (e jT -e kT ) 2 

+2K ija M T (uj) y/ {jicr-jjo){(-j<j-tia) \/ (/fcr-/;r)(e;r-e/cT) 

(8) 

where ti a is the KS eigenvalue of the i-th KS orbital of spin a. fi„ is 1 if the i-th orbital of spin a is occupied; 
otherwise, fi a is 0. uii is the 7-th excitation value sought. Kij a M T {uS) depends on the exchange correlation kernel as 
follows: 



/ 



dx dx'(t>* a (x)<t>* jlT (x) f KXOi<rr (x,x',u) (]) kT {x')<j)i T {x'). 

(9) 



In general, the matrix K(w) depends on the energy, u>, of the solution, and this greatly complicates the solution of 
the general problem. However, the u> dependence is often ignored in what is called the adiabatic approximation. 

The adiabatic LDA kernel can be derived from the second functional derivative of the time-independent LDA 
exchange correlation functional. Hence, 

(LDA ( I \_ S2E Hxc A i n i(] _ 
JHXC ar\ x i x ; u I 



Sn (T (x)6n T (x') 

X6(x - (1 - S„) + f ^fcHx))} + 2f£ (n(^ _ y 
\ d^nyx) n(x) ) 

= f^ x aT {x,x') + f^{x,x'). 

(10) 

The kernel is split into two parts, an Hartree-exchange part (/ffx^r) an( ^ a correlation part (f^f). The former is 
known exactly, the latter must be approximated. The result after the second equal sign is tailored specifically to the 
application in this manuscript assuming an unpolarized ( — system. 

Now, we will explain why we need to have an accurate correlation energy per-particle to third order in the high- 
density limit. The second derivative of the correlation energy density with respect to the density represents the 
heart of the adiabatic TDDFT approximation as it carries all the correlation effects beyond what is modeled in the 
KS orbitals. The correlation kernel is the second functional derivative of Eq. 6 with Eq. 5 plugged in explicitly. 
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In the high-density limit, the first two terms of its second derivative vanish, and only the third and higher terms 
remain. Therefore, in order to model the correlation kernel at all in the high density limit, the third term must be 
included. For this reason, we had to replace the (3,3) forms with the (4,4) Pade given in this paper. In 3D, the 
situation is different as the logarithmic dependence on the density means that the second derivative is divergent and 
non- vanishing; however, any additional finite terms might be neglected, and third order non-logarithmic terms might 
become important in certain common density ranges. 

The ID H2 model at hand has at most 2 bound KS orbitals, so we use the two-state single-pole approximation to 
TDDFT, first presented by Casida in Ref. [22]. It is assumed that the contribution of any finite number of unbound 
orbitals to the response function is negligible due to box normalization. It is possible however, in the limit of extremely 
weakly bound orbitals, that this approximation is no longer valid. The solution of Eq. 8 within the two-state model 
is 

LO S = y/(ei - e )[(ei - Co) + 2(K n + K n )] (11) 

for the singlet excitation and 

lot = yftei - eo)[(ei - e ) + 2(K n - K n )] (12) 

for the triplet excitation. e corresponds to the lowest-energy gerade KS orbital eigenvalue and t\ corresponds to 
the ungerade excited KS orbital eigenvalue. We use eo and t\ here to distinguish these values from the exact non- 
interacting eigenvalues e+ and e_ even-though the KS states are also characterized as gerade and ungerade. K is a 
2x2 matrix in spin given by 

K aT = \{l-5 aT ) dx4>l(x)4>l{x) 

/CO /'OC 
/ dx dx' Mx)Mx)fc D ar(^ X '^)Mx')h(x')- 
OO J — OC 



(13) 



In section V, we will directly compare the exact excited-state spectrum for Eq. 1 and the TDDFT spectrum. 

DFT will likely fail for large separations, a, because the local treatment of correlation does not cancel the exchange 
and Hartree terms. This is the long-ranged self-correlation error, sometimes called the static correlation problem in 
DFT, and it exists even when the interaction is local. The net effect is that two well separated fermions interact in 
the LDA while in a realistic system, the fermions would be entangled but otherwise non-interacting. For ground-state 
DFT, we suggest a simple scheme to model the long-range correlation and to cancel the self-correlation error in the 
stretched case. First, note that in the stretched situation, there is no problem in the polarized case since two fermions 
of the same spin do not interact via the contact interaction. In the unpolarized case when £ = 0, opposite spin 
particles are likely to localize in different regions of space, the interaction energy will be much less than expected if 
they delocalized. In order for this to be captured, the correlation energy must exactly cancel the Hartree-exchangc 
energy. This is achievable if ec(ji) = —An/4 when £ = 0, the negative of the Hartree and exchange energy. The 
solution is then to obtain information from the density that fermions are in the stretched regime and to apply long- 
ranged correlation in this case. A previous attempt at describing this long-ranged correlation used the pair density 
function with some success [23]. Another way this can be achieved is by relying on the KS orbitals and defining the 
dimcnsionless and unitary parameter, 



E oc J<M 4 Qe) 
(E occ . I<M 2 W) 2 



<*) = £<?%%v > ( 14 ) 



where <pi t(T represents the i-th KS spin-orbital, and the sum is over i and a. If the fermions are isolated as in the case 
of a one fermion system, t(x) = 1; otherwise, t(x) < 1. For a two fermion system such as Diracium, t(x) = 1/2, 
and t(x) — for the uniform reference system. Physically, r can be thought of as a measure of the inverse number 
of fermions that are locally relevant, and in this way, transcends the definition given here. Unitary r implies that a 
certain region of space is occupied by only one fermion and so many body effects are trivially unimportant. As r gets 
smaller in magnitude, the importance and nature of many body effects becomes important. 

We restrict ourselves to the case £ = 0. Suppose that we use t(x) to model a local toggling between uniform- 
reference-system-based DFT and the exact long-ranged limit: 



/oo 
dx [f(T(x))e^ DA (n(x),0)-(l-f(T(x)))e H x(n(x),0)] n(x). 
-OO 

(15) 
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/(t) = if t = 1 and otherwise has a value of unity for the two fermion system, but in general, a more complicated 
form is necessary. We call this functional LC for locally cognizant referring to the fact that the approximation is 
similar to the LDA but allows for the recognition and proper treatment of single particle regions. Eq. 15 produces 
the same LDA results for Diracium and the small-spacing double- well problem since t(x) = 1/2 everywhere for these 
systems. For stretched ID H 2 , t(x) = 1 for the exact solution. LDA gives t(x) — 1/2 and the wrong energy. 

To handle the large separation limit, we need to consider details about the Kohn-Sham reference system. In this 
limit, the highest occupied and lowest unoccupied approximate KS orbitals are nearly degenerate, and it is plausible 
that the exact KS potential would result in orbitals that most closely resemble linear combinations of these. If we were 
to express the exact KS orbitals in terms of the approximate ones, we would no longer be in a variational minimum of 
the approximate KS system, and consequently, the orbital energy contribution would rise. Our approximate density 
function must then lower the interaction energy to compensate for this effect. 

This LC functional will give the correct result if the interaction energy is lowered by a greater amount than the 
Aufbau rule raises the energy. The scheme implies that symmetry can be broken. There are two lowest energy KS 
solutions: one with the up fermion localized right and the down localized left and vice verse. The total density for 
each is the same, but the magnetization is inverted. Since both solutions are of equal weight, the physical observable 
are expected to be averages of the equal energy states. Thus, the total magnetization will vanish. The philosophy 
here differs from the wave-function based idea of multiple determinants. 

For example at a = 2, the energy change of elevating the occupied orbitals is AT S w 0.037, and the change of 
turning on the long-ranged part, AE C = E^ DA — E HX = —0.157. These two conspire to lower the total energy. For 
o=l, the AT S = 0.259 dominates over AE C = —0.195, and the local correlation method is valid. Examination of Eq. 
15 shows that LC will over-correlate relative to the LDA but will also increase the kinetic energy. We will see later 
that this performance is required to improve upon the pure LDA. An optimized effective potential scheme is needed 
to apply LC self-consistently, even to this ID system and will be explored in future work. The stretched limit poses 
challenges for TDDFT as well. In the large separation limit, we should have fc a,r(x, x' , u>) = —\5(x — x') (1 — 8 aT ) 
because in the large separation limit E c — —E HX . The form of Eq. 10 is compatible with this limit if ec = — An/4 as 
the exact correlation functional should be in this limit. It would be interesting to known whether the self-interaction 
corrected LDA would correctly give this limit. This too requires an optimized effective potential scheme. 



IV. METHOD OF EXACT SOLUTION 

Here, we present a technique to exactly solve the eigenvalue problem, H^/(xi, x%) = E^(xi,X2) using H in equation 
1. This work extends an idea originally introduced in Rosenthal's work [13] but involves many more challenges than 
the Diracium solution. 

The exact spin-singlet real-space wave-function of two ID fermions is a 2D function, ip{x, y), constrained by particle 
interchange rules to be either symmetric or antisymmetric under the swapping of x\ and x^. The eigenvalue problem 
is difficult to solve using the traditional wave-function based methods of quantum chemistry such as configuration- 
interaction, because the virtual spectrum is mostly unbound. A prohibitively large number of excited configurations 
would be needed to provide an accurate solution. Perturbative approaches suffer similar limitations. A variational 
method could provide a highly accurate solution. However, the solution would be affected by the assumed form of the 
variational wave-function. This form, in principle, limits the accuracy of the many-body solution. In this section, wc 
present an exact numeric solution. By exact numeric, we mean an integral equation that can be solved to arbitrary 
accuracy numerically by increasing and refining the number of integration points. 

The exact solution is found by reducing the Schrodinger equation with the Hamiltonian in Eq. 1 to a set of coupled 
integral equations. To do this, the differential equation is expressed in momentum space in terms of three ID trace 
functions: 

/oo 
dy exp(iky) [ip(a, y) + ip(-a, y)] , (16) 
-oo 

/oo 
dy exp(ifcy) [ip(a, y) - ip{-a, y)] , (17) 
-oo 

and 

/oo 
dy exp(-i%)V(y,y). (18) 
-oo 
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Note that we use y here as the conjugate spatial coordinate to fc to avoid confusion with x\ and x 2 . The fc-space 
wave-function solution can be expressed algebraically in terms of these trace functions: 

$(h,k 2 ) = 
2 

2 (cos(fcia)Gi(fc 2 ) + cos(fc 2 a)Gi(fci) 



k 2 



+ sin(fcia)G 2 (fc 2 ) + sin(fc 2 a)G 2 (fci) - XH (fci + fc 2 )). 

(19) 

The many-body energy eigenvalue is E = —p 2 /2 defining p. 

To obtain a set of ID integral equations, we Fourier transform Eq. 19 into real space and use the result to express 
ip(a,y) and ip(—a,y). Then, we plug ij){a,y) and ip(—a,y) into the right-hand sides of Eqs. 16 and 17. 

For So, the ground-state, we write out the resulting coupled integral equations for Gi and G 2 explicitly, 

2 (, (l + exp-V( afe ) 2 +( a P) 2 )\ _1 



GUak) = - 1 



d(afc') 



2a cos(afc) cos(afc') 
(afc) 2 + (afc') 2 + (ap) 2 



2Aa 2 



7T 

2Aa 2 



{n C c,i,ap,a\{-ak, -ak') + K C c.i,a P ,ax(-ak, ak')) Gi(afc') 



and 



(K cs ,i,ap,ax(-ak,-ak') - K, CStl , aPia \(-ak,ak'))p 2 (ak') 



2 ( (i + exp -V(^) 2 +(^) 2 ) 

G 2 (afc) = — 1 — a 



(20) 



7T 



d(afc') 



2asin(afc) sin(afc') 



(afc) 2 + (afc') 2 + M 



2Aa^ 



7T 

2Aa 2 



{K ss ,i,a P ,a\(-ak, -ak') - n as ^^ a \{-ak, ak')) G 2 (afc') 



(K sc ,i,a P ,ax(-ak, -ak') + K sc ^^ a \{-ak, afc'))Gi(afc') 



(21) 



with 



K cc,a,p,\(k, k ) — 



CM 1 



V"? 2 + 2 p 2 



/ cos[a(fc + q)} \ / cos[a(fc / + g)] \ 
V(fc+g) 2 + fc 2 +p 2 / \(k'+q) 2 + k' 2 +p 2 ) 



+ 



cos[a(fc — q)] \ / cos[a(fc' — g)] 
(k-q) 2 + k 2 +p 2 ) \{k'-q) 2 + k' 2 +p 2 



(22) 

We have placed pre- factors of a in locations that are convenient for numerical reasons. A similar set of equations 
can be written down for the Si state. S 2 satisfies the above set with a different p. K cc , a ,p,A(fci, k 2 ) can be evaluated 
explicitly using complex analysis. The notation cc stands for the trigonometric functions that are contained in k. cc 
is for cos cos, cs for cos sin, and so on. The square-root introduces a branch cut affecting the contour integration that 
is best handled through Gauss-Legendre numerical integration. 

So, we arrive at a set of coupled ID integral equations in two functions Gi(fc) and G 2 (fc). Converting the integrals 
to Gauss-Legendre sums allows us to express the coupled integral equations as a matrix problem: 



Gi 
G 2 



Mn M 2 i 
Mi 2 M 22 



Gi 
G 2 



(23) 
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State 


nonint 


occ. 
1 


orb. 

2 


exchange symmetry 
space spin 


two-body 
energy 


So 




+ 


+ 


Symm. 


Asymm. 


2e+ 


Si 




+ 




Symm. 


Asymm. 


e++e- 


s 2 








Symm. 


Asymm. 


2e- 






+ 




Asymm. 


Symm. 





TABLE I: The non-interacting many-particle bound-states of the Hamiltonian Eq. 1. + stands for the gerade single particle 
orbital with eigenvalue e+, and — stands for the ungerade single particle orbital with eigenvalue e_. 



where n(p) is an eigenvalue that equals unity when the appropriate p, that satisfies the original eigenvalue problem is 
inputted. 

In expressing the integral equation in matrix form, the integrals are discrctizcd onto a set of points. There is no 
unique way to do this. We chose a Gauss-Legendre inspired integration scheme. However, the implementation is not 
straightforward as the integrals take on the form, 

/°° trier 2 
J k WTwwrm m - <24) 

where the function trig 2 fc can be cos 2 k, sin 2 k, or sin k cos k and £(fc) is a smooth non-oscillatory function of k. To 
obtain accurate numerical quadratures, we break up the domain into intervals of tt/2 and integrate each separately 
using a suitable method. For convergence to better than 10 nano-Hart. in energy and 10 -8 in at a = 1, we need 
800 grid points for the primary integration region and 200 additional tail points to model a portion of the asymptotic 
tail. The complicated scheme is highly accurate and has been verified by producing the exactly known non-interacting 
results to 9 significant figures for 1000 k-points. 

Two tests verify the accuracy of the exact solution. First, the exact solution is shown to approach the known 
united-atom limit. The idea is that this exact solution should give results approaching Diracium's with an appropriate 
combined potential strength. For two wells with Z = 1 each, we get Diracium with Z = 2. Figure 1 shows the energy 
for ID H 2 for various a and methods, and the approach can be seen. The approximate methods will be discussed in 
section V. The far left limit on the plot is Diracium when A = 1 and Z = 2, _E 2 = —3.155. For crushed iJ 2 , we find 
E = —3.023 when a = 0.01 and E = —3.087 when a = 0.005 extrapolating to E 2 = —3.152 in excellent agreement 
with the united-atom limit. 

As a second test, the exact numeric solutions are shown to satisfy the generalized virial relationship given by Eq. 
2 to within 10 mHart. The accuracy is modest due to the numeric challenge of solving the required triple integrals 
with the limited sampling of k points. For Z = 1 and A = 1, the left-hand side of Eq. 2 is —318 ^Hart. while the 
right-hand side is —317 /xHart. in less than perfect but still acceptable agreement. 



V. COMPARISON OF TDDFT AND HUBBARD MODELS TO THE EXACT RESULTS 

In this section, we compare adiabatic TDDFT and the Hubbard model to the exact results. To start, we make 
some comments about the parameter regimes chosen. This model offers a rich spectrum of phenomena; however, for 
the sake of brevity, we have restricted ourself to a rather arbitrarily chosen parameter ranges. As mentioned, we set 
Z = 1 in this section for convenience, but results for arbitrary Z can be related via scaling. 

We analyze the energy spectrum of ID H 2 versus A for two values of a. The first case is the double- well potential 
with 2 bound- fermions and separation a = 1. In this case, we expect DFT to provide an accurate description of the 
spectrum since the fermions are both localized in the area of the double well. The Hubbard ground-state is expected 
to be too low by at least 68 mHart, the amount that the kinetic energy is misrepresented in the non- interacting case. 
For the second case a = 2, the system is in the stretched H 2 limit. The LDA is known to be unreliable in this limit 
because its failure to capture the localization of fermions to opposite sites without symmetry breaking. On the other 
hand, the Hubbard model is designed to work well in the stretched case. For example, the non-interacting Hubbard 
error at a = 2 is less than 3 mHart. for the ground-state, So- For now, our primary focus is on comparing stable 
results, and thus, we restrict our analysis to A < 2, safely within the bound regime for the ground-state, So- 

Our labeling scheme for the quantum states is motivated by the noninteracting many-particle spectrum. The 
non-interacting multiple particle spectrum can be constructed from the single particle states. For the 2-particlc 
solution, products of the single particle orbitals must be properly symmetrized. The entire bound spectrum of the 
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FIG. 1: (Color online) Normalized ground-state energies for the ID analog of H2 at separation a and interaction strength, 
A = 1, within several approximations. The results are normalized by the noninteracting result. In the asymptotic limit, the 
ratio should be unity. The solid line is the LDA result, the dotted line is the LDA LC result, the medium-dashed line is the 
Hubbard result, the long-dashed line with diamonds is the exact result, and the arrow indicates the united-atom Diracium 
limit. 



non-interacting 2 particle states is given in Table I. Note that the non-standard labeling of the states is used because 
angular momentum needed for the traditional labeling scheme is not well-defined in ID. There are 3 spin-singlet states 
and 3 spin-triplet states. The spin-singlet states have symmetric wave-functions under particle interchange. There 
are 2 single-particle orbitals, gerade and ungerade ones, and 3 unique symmetric products can be made of these. 
The triplet state is triply degenerate in the spin manifold but has only one spatial contribution, an anti-symmetrized 
product of the gerade and ungerade orbitals. Antisymmetric products of like noninteracting orbitals vanish. It is 
pointed out that the singlet-state comprised of a product of 2 ungerade orbitals represents a double excitation from 
the ground-state. The states are labeled So, Si, and S2 in the singlet manifold, and T\ in the triply degenerate 
triplet manifold. A transition form So — ► Si is the first singlet excitation, and a transition form So — > 2i is the first 
triplet excitation. Both are calculable in first-order linear response theory of DFT. A transition form So — > S% is a 
double excitation, proportional to the intensity of light squared, and is formally beyond first order response theory. 
This fact can be realized by considering the noninteracting response function and noting that no double poles exist. 
Nevertheless, the bare Kohn-Sham (KS) double excitations might be considered a good first approximation to the 
double excitation due to the linear response terms vanishing and a better approximation than the KS single excitations 
were to their counterparts. 

The interacting two-fermion spectrum is expected to have a one-one correspondence with the noninteracting one 
excepting the possible disappearance of the highest energy states into the continuum as the interaction, A, is increased. 
Take for example, the ID analog of Helium, Diracium, at Z = 1 . A single 5- well typically can bind two fermions of 
opposite spin in ID, but above A c «t. = 2.6673532258, the two-particle state merges with the continuum, and only 
one fermion can be bound [13, 14]. For two- fermions in a double- well, the critical interaction strength, X cr it.H 2 { a i S), 
depends on the well-spacing, a, and state, S. The critical interaction strength, \ C rit.B 2 ( a i is likely larger than its 
corresponding value for a single well due to the stabilization effects of the hybridized orbitals or, in Hubbard parlance, 
the energetic favor-ability of hopping. Most likely, each excited state in the double well has a different critical value 
as the higher energy states are likely to vanish at smaller A than the ground-state. However, the triplet state is 
interaction independent because like spins do not experience the contact interaction. 

In Fig. 1, the total energies are normalized by the non-interacting results for A = 1. In the limit of large a, the 
ratio should become unity as the well-separated fermions will not interact but will reside on different sites. The 
exact solution is globally spin unpolarized but locally acquires a nontrivial spin dependence. This is the statement 
that the two fermions will occupy different sites and break symmetry. The exact Sq result here is the fully non- 
constrained solution of the ID H2 Hamiltonian. The Hubbard model Sq reproduces the long range charge separated 
limit accurately but fails to describe the crushed limit due to improper scaling of the hoping term in this limit. This 
is easily seem through the non-interacting scaling behavior of the hoping term. The LDA/EXX So does approach the 
crushed atom limit as is expected from earlier work on the Diracium system. However, the LDA alone fails to capture 
long-ranged charge separation and fails to capture about 25% of the total energy in the stretched limit. A perturbative 
application of the LDA LC So greatly improves this approach since it allows for opposite spins to become separated. 
In the intermediate range, a ~ 1 — 2, the DFT methods and Hubbard model bracket the exact result indicating that 
both proper scaling and ability to localize particles are vital in this range. 

In Fig. 2, we see the energy spectrum plotted for ID H2 with a = 1 at various interaction strengths, A. The exact 
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FIG. 2: (Color online) DFT and exact energy spectrum in atomic units for the ID analog of H2 for various A with inter-atomic 
separation, a = 1. These results are generated through the exact solution of the quantum many-body problem and the LDA 
version of TDDFT. So and Si are the ground-state and first singlet excited state, and Ti is the first triplet excited state. The 
solid line is the LDA So result, the gray varyingly dashed line is a restricted version of the exact result, the long-dashed line 
with diamonds is the exact result for So, the dotted line is the LDA Si result, the medium dashed line is LDA Ti result, and 
the alternating short-dashed line is the exact Si result. 




FIG. 3: (Color online) Hubbard and exact energy spectrum in atomic units for the ID analog of H2 for various A with inter- 
atomic separation, a = 1. The long-dashed line with diamonds is the exact result for So, the alternating short-dashed line is 
the exact Si result, the long-dashed line is the Hubbard So, the alternating short-long line is the Hubbard Si result, and the 
alternating medium-dashed line is the Hubbard Ti result. 



result does not level out with increasing interaction strength. The system probably ionizes at some critical interaction 
strength as was the case for Diracium. Exact ground-state (Exact So) and restricted exact ground-state (Res. Exact 
Sq) results are presented. In the later, the solution is forced to be expressed purely in terms of G\(k) with the 
antisymmetric G2{k) forced to vanish (See Eq. 20 and 21) . In wave-function theory, this corresponds to unbroken 
symmetry. The result is a higher energy solution than the true broken symmetry ground-state. The difference for 
weak interactions is negligible but at larger interactions the restriction causes an energetic error exceeding 10%, and a 
qualitative prediction of the cross over between the triplet excited-state and the ground-state. The exact result does 
not cross over. This finding highlights the importance of symmetry breaking in describing static correlation in the 
design of density functionals. 

LDA/EXX So and the restricted exact So results are in excellent agreement until A = 1. Beyond this point the two 
deviate slightly but increasingly. This is most likely due to the difficulties of approaching the ionization threshold 
where the density is less localized. The S\ states agree up to A = 0.5. The LDA S\ fails to display the leveling off of 
the exact S\ result. The DFT triplet (LDA/EXX Ti) result like the Hubbard (shown in Fig. 3) and exact results is 
constant. However, at A = 1.2, the triplet state becomes unreliable. This is due to the cross over between the DFT eo 
energy and the t\ orbital energy. The single pole approximation becomes numerically unstable and thus unreliable. 

For the ground-state, KS theory performs quite outstandingly reproducing 98 % of the restricted ground-state 
energy at A = 0.5. In the strongly-interacting regime at A = 2, DFT still gives a result within 80 % of the exact value. 
EXX on the other hand would be more inaccurate giving about 60% of the total energy. The EXX method is related 
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FIG. 4: (Color online) Hubbard and exact energy spectrum in atomic units for the ID analog of H2 for various A with inter- 
atomic separation, a — 2. The long-dashed line with diamonds is the exact result for So, the alternating short-dashed line is 
the exact Si result., the alternating medium-dashed line is the Hubbard So, and finally, the alternating long-dashed line is the 
Hubbard Si result. 




FIG. 5: (Color online) DFT and exact energy spectrum in atomic units for the ID analog of H2 for various A with inter-atomic 
separation, a — 2. The long-dashed line with diamonds is the exact result for So, the alternating short-dashed line is the exact 
Si result, the dotted line is LDA LC So result, the medium-dashed line is the DFT Ti result, and the long-dashed line is the 
DFT Si result. 



to the Hartree-only theory and the nonlinear Schrodinger equation that is popular in many treatments of ID systems. 
We see here an example that non-linear Schrodinger approach does not work accurately or reliably for moderate 
to strong interaction strengths. Although the DFT results perform exceptionally well up to moderate interaction 
strengths, the results are under-correlated for stronger interactions, A > 0.5. This is because partial localization 
of the fermions to opposite ends of the system is not accounted for. The exact restricted results correspond to a 
solution of the exact integral equations with G2(k) forced to vanish. This is valid for the non- interacting solution 
because in that case the eigenstates are single Slater determinants of non-interacting single-particle solutions. For the 
interacting state, the restricted result does not represent a true eigenstate of the original Hamiltonian and according 
to the variaitonal principle, has a higher energy than the exact ground-state. Notwithstanding, the restricted energies 
agree quite well with the pure LDA values indicating that LDA correlation is adequate to describe the system if there 
were no localization. 

TDDFT allows us to find the spectrum of excited-states. For the Si state, we show only up to A = 1 in Fig. 2. 
At larger interaction strengths, this state tends to be unstable and decays into an unbound Fermion and one bound 
Fermion. As a measure of this accuracy for TDDFT in the two state model, we work backward. The triplet state 
does not experience the interaction so the triplet excitation is known exactly without using the TDDFT formalism. 
The triplet is unaffected by the interaction and should be a straight horizontal line with respect to A. Thus, TDDFT 
should reproduce this line if the kernel and orbitals are both accurate. This is what is seen up to about A ~ 1.2. At 
this point the triplet energy and the singlet restricted ground-state energy are close hinting at a level crossing. The 
occurrence of a cross-over differs from the Hubbard model where no level crossing occurs. The total energy excited 
state gap for the singlet is much larger than for the triplet gap. It is interesting to note that the singlet gap at 
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A = 1 is only about one quarter larger than at A = 0. Adiabatic effects are important, but the ratio of the singlet to 
triplet gaps grows more than the energy of the singlet state. The actual non-adiabaticity must be a very complicated 
functional of the excited-state energies to capture this behavior. For the triplet to be well reproduced implies that the 
potential, orbitals, and kernel (at this energy range) are accurate. But the prediction of the first singlet excitation is 
not as accurate. Since the orbitals and potential are the same for the singlet, the approximation of the kernel must 
not be as accurate in the calculation of the singlet. 

The Hubbard So is typically over-correlated as shown in Fig. 3. This is, in part, due to overemphasis of the kinetic 
energy at this smaller spacing. But the model is also qualitatively wrong for larger interactions as the energy goes 
to a constant while the exact result continues to grow with the interaction strength until the system is eventually 
ionized. The over-correlation is due to improper scaling towards the crushed H2 limit. Additionally, the two-site 
Hubbard Hamiltonian does not allow for unbound-states and is consequently incapable of describing ionization. In 
Fig. 3, the ground-state (So) Hubbard curve resembles the exact curve except for a offset of about 0.1 Hartree. This 
lower energy is a result of the improper scaling of the kinetic energy hopping term in the united-atom limit. The 
Hubbard Si state suffers from a similar offset problem, however, the Si state does not demonstrate the sharp leveling 
off of the exact result at A = 0.6. The higher energy Hubbard result can be rationalized by realizing that for large U, 
there will be significant projection of the localized solutions onto the other site thus increasing the effective overlap. 
This increased derealization is not described by the Hubbard model. The T\ Hubbard curve is flat by construction 
since it does not depend on U, the inter-particle interaction. 

For the stretched case, we expect the Hubbard model to be essentially exact because the hopping term and single 
site repulsion terms are accurate. This is in fact what in seen in Figure 4. The remarkable agreement is not surprising 
as the single orbital overlaps decay exponentially, and the relevant hopping parameters are small. The 7\ Hubbard 
curve is again level as explained in the previous plot. For larger A, the Hubbard result for Si is dangerously close to 
ionization w —0.5. It was seen that local correlations are inadequate to properly describe ionization. 

In Figure 5, the comparison between DFT and exact is less satisfactory. The ground-state results deviate quite 
substantially for interactions just larger than A = 0.1. The ground-state energy is in error due to the failure of LDA 
to account for the localization of fermions to opposite sites. The LDA LC, applied perturbatively, does provide some 
of this information and thus drastically improves the agreement with the exact 5*0 result. A self-consistent application 
of this functional would involve an optimized effective potential algorithm that is beyond the scope of this work. 
However, self-consistency is likely to improve the accuracy. It would be interesting to find out how a kernel based on 
this improved functional would perform. The perturbative approach only approximates the stretched H2 limit while 
a fully self-consistent approach should exactly match at large A. LDA/EXX Ti only agrees up to to A = 0.1. The 
orbitals are not faithful representations of the exact KS orbitals since the LDA ground-state is inaccurate. Curiously, 
the LDA/EXX Si and exact Si agree exceptionally well up to A = 1/2. This is because the local kernel cancels the 
self-interaction correlation error in the ground-state calculation. No excited state results are reported for LDA LC 
TDDFT because the method had been applied perturbatively, and the self-consistent wave-functions and kernel are 
not available in the analysis. 

VI. CONCLUSION 

In this paper, we have explored how two standard viewpoints of condensed matter physics describe an interacting ID 
many-particle system. One method, LDA DFT, provides accurate energies within about 0.1 mHartree for intermediate 
interaction strengths (A < 0.5) and distances (a < 2) but fails appreciably at large well spacings. Analysis of the exact 
result shows that this limitation is due to symmetries induced by the short-sightedness of the restricted KS scheme 
within the LDA. When the restricted symmetry is enforced on the exact solution, the result lies much more closely 
to the LDA indicating that the lack of localization is the key deficiency in the restricted LDA KS treatment. To 
overcome this challenge, a local parameter, t(x), is introduced that describes the effective local number of fermions 
and is readily implementablc in existing electronic structure codes. This local measure allows the use of two reference 
systems in the construction of a density functional: the uniform reference system and the single particle system. The 
introduced functional, when applied perturbatively, is shown to better reproduce the energy curve of ID H2 versus 
well spacing. A self-consistent application will require an optimized effective potential approach beyond the scope of 
this work but is likely to improve the agreement. 

On the other hand, the two-site Hubbard model provides a qualitatively accurate description of the ground-state 
across a wide range of parameters describing both united-atom and separate-atom limits, but it fails in its quantitative 
predictions and has questionable scaling characteristics. This is not surprising as the model is limited by design. The 
2-site Hubbard model does not include the continuum and will therefore fail to describe ionization and scattering. 

The excited-state results follow a similar pattern as the ground-state results. DFT is accurate for small well 
spacings, and the Hubbard model is more reliable for larger spacings. In TDDFT, the approximation of the higher 
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FIG. 6: Diagrammatic representation of the third order contributions to the interaction energy. Spin labels are omitted since 
the two loops in each diagram must have opposite spins. Up arrows represent particles, and down arrows represent holes. The 
final diagram contributes twice because a horizontal rotation produces a new diagram. 



r 

► 




q 



FIG. 7: Momentum labels for the third order contributions to the interaction energy, q and r are momentum transfers, fci and 
can label particle or hole momenta depending on the diagram. 



energy singlet fails for weaker interactions than the lower energy triplet. For the triplet to be well reproduced implies 
that the potential, orbitals, and kernel (at this energy range) are accurate. Since the first two are the same for the 
singlet, the approximation of the kernel must not be as accurate in the calculation of the singlet. A strongly non- 
adiabatic kernel would explain why TDDFT predicts the lower energy triplet but not the higher energy singlet. Thus, 
non-adiabaticity outweighs ultra-nonlocality problems. Perhaps, this is due to the locality of exchange for contact 
interactions and could be quite different from what occurs with long-ranged 3D interactions. Curiously, the Hubbard 
first excited-state can also prove drastically wrong in cases when the fermions are strongly interacting, large U, as 
the Hubbard treatment forces same site localization while the exact system will have significant delocalization. This 
result has significant implications for the reliability of LDA+U results with large U. 

The realm of ID contact-interacting fermions provides an interesting opportunity to compare the Hubbard and DFT 
models in detail. Results found here have provided insight into models of 3D Coulomb interacting systems where 
many of the underlying quantum many-body effects are obfuscated by the long-ranged nature of the interaction. 



APPENDIX A 



In this appendix, we find the third-order A term in the high density limit of the correlation energy per particle 
for Deltium, the one-dimension uniform fermion system, using the Goldstone diagrammatic approach to perturbation 
theory in momentum space [27]. 

The Fourier transform of the interaction potential is V{q) = ^ J^ L dx 5(x)e lqx = ^ where L is the arbitrary length 
between the boundaries confining the system. Like spin fermions do not interact via the 5-function, this means that 
only vertices that connect opposite spins enter into the diagrammatic series. This is a tremendous simplification as 
many diagrams vanish. A further simplification is that the interaction is independent of the momentum transfer, q. 

To third order, three different diagrams contribute as seen in Fig. 6. The momenta in the each diagrams are labeled 
according to the arrows in the two-bubble diagram shown in Figure 7. From the standard rules of perturbation theory, 
the third order term can be written as a sum of multi-dimensional integrals, 
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with |fci + <?| > fcir, |fci + r| > 0, |fc 2 — q\ > k F , and \k 2 — r| > for the first term, with \ki + q\ < k F , \ki + r\ < 0, 
1^2 — q\ < k F , and \k 2 — r\ < for the second term, and with \ki +q\ < k Fl \k\ + r\ < 0, \k 2 — q\ > k Fl and \k 2 — r\ > 
for the final term. The notation d\k\ stands for the sum of two integrals, fj^ dk + dk. k\ and k 2 are particle 
(or hole) momenta, and q and r are the momentum transfers. n s is the number of spin species, in this case 2. The 
limits of integration and constraint inequalities ensure that particle have less momentum than the Fermi-momentum, 
and holes have higher momentum than the Fermi momentum. There is a symmetry factor of 1/2 associated with each 
diagram. The third diagram contributes twice because a horizontal rotation produces a new diagram with the same 
numerical value. This third diagram contributes a negative value because a odd number of vertices connect particle to 
holes. To solve Eq. (A-2) exactly, we re-scale as follows: q = k F v, r = k F x, k\ = k F y, and k 2 = k F z. The correlation 
energy per particle in third order in A is found to be 
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vx(v + y — 


z)(x + y- 


z) 




1 




vx(v + y — 


z)(x-y + 


z) 




1 




vx(v + y — 


z)(x + y- 


z) 




1 




vx(v + y — 


z)(x + y- 


z) 



8tt 2 In 2 -36 C(3), 



1 



00 pv pl+x pl+x 

A I dv I dx I dy I dz- 

h Jv-2 Jv-i Jv-i vx(y + z-v)(y + z~x) 

r 2 r v rl+x rl+x 



pZ pv pl+x pl+x 

+4 dv dx dy dz 
Jo Jo Ji Ji 



vx(y + z — v)(y + z — x) 
= -87r 2 ln2 + 48 C(3), 



and l c 



p2 pv pl+x pi ^ 

4 dv dx dy dz — j— — 2 

Jo Jo Ji Jl-x vx{y-\-z) 



1 



/•4 pv pl-\-x pi 

H dv L dx L dy L dz v X (y+zr 

r4 rv rl-\-x f>l -| 

+41 dv j dx j dy j dz —j- ■ 



v—2 Jv~l Jl — x 



vx(y + z) 2 



(A-2) 



(A-3) 



(A-4) 
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oo pv ^l+a; pi ^ 



+4/ dv dx dy dz 

Ji Jv-2 Jv-i J -i vx(y + z) 

= 2C(3). 

(A-5) 

The author suspects that an even more tedious calculation reveals that the correlation energy per particle to fourth 
order in A is 

C W = = -0.00093?4. (A-6) 

C 47r 6 n 2 n 2 

This value agrees well with the numeric Bethe-Ansatz result of -0.00094 X 4 /n 2 . The full derivation will be given in a 
later work. 

APPENDIX B 

In the (4,4) Pade parameterization of e xc ("-, C)i we have used the result that there exists an astonishing non-linear 
analytic relationship between the given expansion limits and the desired Pade parameters. A (4,4) Pade of the form 

„. . Ax* + Bx 3 + Cx 2 

F{x) = flWT^ (C " 1} 

has parameters A, B, C, D, E,a,nd F chosen to satisfy known limits. If three terms in both the large and small x 
limits are known as given below, we can fit these parameters exactly. Specifically, if for large x, 

F(x) =gix + 92+93,/x + ... (C-2) 

and if, for small x, 

F(x) = g 4 x 2 + g 5 x 3 + g 6 x 4 + .... (C-3) 

the parameters A-F can be determined explicitly. 

Let us introduce the following nonlinear Ansatz for the parameters. 

A = .91 (.94 + 2.91.95.94 + 323654 - ff2ffl + gfae) 



9l + 2 ff2S'4ffl + 535551 + 5354 - 5255 



B = 555? + ghi 



(CM) 



+ 



5i + 2525451 + 535551 + 53.94 - 5255 
(525455 + 53 (fff ~ 545e))5i+52(5l + 525654 ~ 52.95) 
5i + 2525451 + 535551 + 5354 ~ 5255 



(C-5) 

C = , 94 , (C-6) 
D = A/ 9l , (C-7) 



and 



E = 55g| + gffll - .92.9651 + 5351 ~ 525455 ~ 535456 

5i + 2525451 + 535551 + 5354 - 5255 



F = gfgi ~ (5255 + 5356)51 + 5254 ~ 535455 + .9256 

5i + 2525451 + 535551 + 53.94 ^ 5255 



(C-8) 



(C-9) 



Direct substitution of Eqs. C-4 to C-9 into the Eq. C-l gives the desired high and low density expansions, Eqs. 
C-2 and C-3. 
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APPENDIX C 



The exact spin-dependence of the exchange correlation energy was alluded to in the body of the paper. Here, we 
present the exact result for the high-density limit. 

For completeness, recall that the spin-dependent Hartree and exchange energy per particle is 

eHx{ni,n L ) = An T n x /(n T + 

In the high-density limit, correlation energy contributes to second order in A. This contribution is described by 
the two-bubble diagram shown second in Figure 7 of Ref. [4] . From the standard rules of perturbation theory, the 
diagram can be expressed as an integral, 

,(2) 



Ne%>(n h n l ) = 
L 2 8tt 3 J_ 00 d<1 J_ kp J_ k ^ q(q + h- fc 2 ) ' 



with kp] = kp[ — n.j_/2, \k\ + q\ > kp^ and \k 2 — q\ > kpi, k\ and k 2 are particle momenta, and q is the 

momentum transfer. Once again the symmetry factor of 1/2 is canceled by the two possible spin configurations. To 
solve Eq. (B-l) exactly, we define two quantities kp and s according to the following: kp^ = kp(l — s) and kp[ = kps. 
Then, we re-scale the coordinates as follows: q = kpx, k\ — kpy, and k 2 — kpz. Note that s = n^/n and £ = 1 — 2s. 
After some algebra, we find the correlation energy per particle: 



#~£(s)!" I <0~5'<0 

4tt 2 

2 

(y-(l-C)HDiLog(l-C) - (l+C)»DiLog(l+C)) 
using the quadrature result below and replacing s by 1/2(1 - (), 

T(s) = / dx dy dz — — ■ 

J2-2s J-s J s-i x{x + y-z) 

p2s f>8 rX+S — l ^ 

+ dx dy dz — — ■ 

Jo Js-x Js-i x(x + y-z) 

2 — 2.S pS rX-\-8 — l ^ 

dx dy dz — — ■ 

2s J-s J s-i x{x + y-z) 

4 ( y - (2 - 2s)5iDiLog(2 - 2s) - 2s5iDiLog(2s)^ . 

The approximation, /(C) w (1 — £ 2 ), is only true for the extreme values of £ with relative errors of up to about 
33.333%. An expansion of ec about £ = for example reveals logarithmic dependencies. In the small (, nearly 
unpolarized limit, 

/(C) = i - 4c 2 + 4> C 2 log C + o(( 3 ). (B-i) 

Perhaps this interesting behavior has implications for the effects of correlation on ID spin density waves, and this will 
be explored further in latter work. Likewise, expansion about £ = 1 shows 

/(C) - 4^DiLog(2C - 2) + 0((C - l) 3 ). (B-2) 
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M 



0.0 0.61478 1.0 oo oo 

0.00068732 0.61459 0.999 42.2794 4.77927 

0.626582 0.44646 0.9 1.19376 1.31123 

1.182951 0.31075 0.8 0.72483 0.93000 

1.720499 0.19480 0.7 0.47586 0.67034 

2.294671 0.09239 0.6 0.28378 0.12800 

2.63967 0.04511 0.55 0.18487 0.29433 

3.04007 0.00886 0.51 0.07634 0.12800 



TABLE II: Sample exact-exchange DFT results through the solution of Eq. D-l for a=l and Z=l. 



APPENDIX D 



The problem in Eq. 1 can be solved analytically within the restricted exact-exchange density functional approach. 
For the ground-state, the solution is spin-unpolarized. The relevant KS equation for one spin-wave-function is 

1 2 
--\7 2 (P(x) + \\cj}{x)\ 2 (j}{x) -Z S{xi ± a)4>(x) = -ecf>(x). (D-l) 

" i=l,± 

This is the non-linear Schrodinger equation in a double-well potential. For simplicity, we only present results for 
o=l. The solution, vanishing at a distance, can be shown analytically to be 

Hx) = [ V / & MJacobiNC (\/5S x,m)\x\<l (D _ 2) 
\ M csch (\/2e(|x| - 1) +x ) |a;| > 1 

where JacobiNC is a Jacobi Elliptic function. The double- well potential forces cusps at the ±1. The constraint can 
be expressed as a transcendental equation. The other constraint is that the wave-function must normalized to unity. 
These two are solved numerically simultaneously. We present some representative numeric solutions in Table II. These 
values and others were used in the paper to validate the LDA code. 
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